# Replication script for the simulations in Leifeld & Cranmer (2018)
# Last change: 2018-01-07


# ==============================================================================
# Preparatory steps
# ==============================================================================

library("statnet")
library("xergm")
library("RSiena")
library("NetSim")
library("ggplot2")
library("gridExtra")

sessionInfo()
# R version 3.3.0 (2016-05-03)
# Platform: x86_64-pc-linux-gnu (64-bit)
# Running under: Ubuntu 16.04.3 LTS
# 
# locale:
#  [1] LC_CTYPE=en_GB.UTF-8       LC_NUMERIC=C              
#  [3] LC_TIME=en_GB.UTF-8        LC_COLLATE=en_GB.UTF-8    
#  [5] LC_MONETARY=en_GB.UTF-8    LC_MESSAGES=en_GB.UTF-8   
#  [7] LC_PAPER=en_GB.UTF-8       LC_NAME=C                 
#  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
# [11] LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C       
# 
# attached base packages:
# [1] stats     graphics  grDevices utils     datasets  methods   base     
# 
# other attached packages:
#  [1] gridExtra_2.0.0      NetSim_0.9           Rcpp_0.12.11        
#  [4] RSiena_1.2-3         xergm_1.8.2          GERGM_0.11.2        
#  [7] rem_1.1.2            tnam_1.6.5           btergm_1.9.0        
# [10] ggplot2_2.0.0        xergm.common_1.7.7   statnet_2016.4      
# [13] sna_2.4              ergm.count_3.2.0     tergm_3.4.0         
# [16] networkDynamic_0.9.0 ergm_3.8.0           network_1.13.0      
# [19] statnet.common_4.0.0
# 
# loaded via a namespace (and not attached):
#  [1] mvtnorm_1.0-5       lattice_0.20-33     muhaz_1.2.6        
#  [4] gtools_3.5.0        digest_0.6.9        plyr_1.8.3         
#  [7] stats4_3.3.0        coda_0.18-1         gplots_2.17.0      
# [10] minqa_1.2.4         gdata_2.17.0        vegan_2.3-1        
# [13] nloptr_1.0.4        texreg_1.36.23      Matrix_1.2-6       
# [16] labeling_0.3        splines_3.3.0       lme4_1.1-10        
# [19] stringr_1.2.0       igraph_1.0.1        munsell_0.4.2      
# [22] speedglm_0.3-1      mgcv_1.8-12         tcltk_3.3.0        
# [25] lpSolve_5.6.13      quadprog_1.5-5      permute_0.8-4      
# [28] MASS_7.3-45         bitops_1.0-6        grid_3.3.0         
# [31] nlme_3.1-128        gtable_0.1.2        magrittr_1.5       
# [34] scales_0.4.1        RcppParallel_4.3.20 KernSmooth_2.23-15 
# [37] stringi_1.1.5       reshape2_1.4.1      ROCR_1.0-7         
# [40] robustbase_0.92-5   boot_1.3-17         trust_0.1-7        
# [43] deSolve_1.12        RColorBrewer_1.1-2  tools_3.3.0        
# [46] mstate_0.2.8        DEoptimR_1.0-4      flexsurv_0.7       
# [49] parallel_3.3.0      survival_2.39-4     colorspace_1.2-6   
# [52] cluster_2.0.4       caTools_1.17.1

set.seed(12345)

numnodes <- 20  # how many nodes should the networks consist of?
timesteps <- 5  # how many consecutive networks?
minval <- -3    # minimum parameter value
maxval <- 3     # maximum parameter value
n <- 100        # number of replications of the comparison experiment
nsim <- 10      # number of simulations for the GOF assessment


# ==============================================================================
# Simulate network sequences from a TERGM and a SAOM process, respectively
# ==============================================================================

# simulation function for TERGM sequences
simul_tergm <- function(density, reciprocity, transtrip, memory, start = NULL, 
    nodes = numnodes, steps = timesteps) {
  param <- c(density, reciprocity, transtrip, memory)
  l <- list()
  m <- list()
  m[[1]] <- NULL
  if (is.null(start)) {
    l[[1]] <- simulate.formula(network(nodes) ~ edges + mutual + 
        ttriple, coef = param[1:3], nsim = 1,  
        control = control.simulate.formula(MCMC.interval = 10000, 
        MCMC.burnin = 10000))
  } else {
    l[[1]] <- network(start)
  }
  for (j in 2:(steps + 1)) {
    m[[j]] <- as.matrix(l[[j - 1]])
    m[[j]][m[[j]] == 0] <- -1
    l[[j]] <- simulate.formula(network(nodes) ~ edges + mutual + 
        ttriple + edgecov(m[[j]]), 
        nsim = 1, coef = param, control = control.simulate.formula(
        MCMC.interval = 10000, MCMC.burnin = 10000))
  }
  obj <- list()
  obj$networks <- l
  obj$parameters <- param
  obj$memorycov <- m
  return(obj)
}

# simulation function for SAOM sequences
simul_saom <- function(density, reciprocity, transtrip, rate, start = NULL, 
    nodes = numnodes, steps = timesteps) {
  param <- c(density, reciprocity, transtrip, rate)
  if (is.null(start)) {
    nw <- matrix(rep(0, nodes * nodes), nrow = nodes)
    diag(nw) <- 0
  } else {
    nw <- start
  }
  nwlist <- list(nw)
  for (i in 1:(steps + 1)) {
    processState <- create_process_state()
    network <- create_network(nw)
    processState <- add_network(processState, network, name = "nw")
    networkIndex <- get_network_index(processState, name = "nw")
    processState <- add_global_attribute(processState, value = 1, 
        name = "timer")
    timerIndex <- get_global_attribute_index(processState, name = "timer")
    modelManager <- create_model_manager()
    time.model <- create_poisson_model(param = 1)
    modelManager <- add_time_model(modelManager, time.model)
    timeUpdater <- create_timer_updater(timerIndex)
    modelManager <- add_time_updater(modelManager, timeUpdater)
    
    effectContainer <- create_effect_container()
    effectContainer <- add_to_effect_container(effectContainer,
        create_effect("density", networkIndex), param[1])
    effectContainer <- add_to_effect_container(effectContainer, 
        create_effect("recip", networkIndex), param[2])
    effectContainer <- add_to_effect_container(effectContainer,
        create_effect("transTrip", networkIndex), param[3])
    
    for (i in 0:(numnodes - 1)) {
      poissonParameter <- param[4]  # NetSim JSS vignette page 13
      poissonModel <- create_poisson_model(poissonParameter)
      changeModel <- create_multinomial_choice_network_change_model(i, 
          networkIndex, effectContainer)
      tieSwapUpdater <- create_tie_swap_updater(networkIndex)
      modelManager <<- add_time_model(modelManager, poissonModel)
      modelManager <<- add_change_model(modelManager, poissonModel, 
          changeModel)
      modelManager <<- add_updater(modelManager, changeModel, tieSwapUpdater)
    }
    
    timeSpan <- 1
    simulator <- create_simulator(processState, modelManager, timeSpan, 
        verbose = FALSE, debug = FALSE)
    test <- simulate(simulator)
    nw <- as.matrix(network)
    nwlist[[length(nwlist) + 1]] <- nw
  }
  obj <- list()
  nwlist <- lapply(nwlist, function(x) network(x, bipartite = FALSE, 
      directed = TRUE))
  obj$networks <- nwlist[-1]
  obj$parameters <- param
  return(obj)
}

# simulate sequences of networks based on TERGM and SAOM process
sim_tergm <- list()
sim_tergm_pos <- list()
sim_saom <- list()
failures <- 0
for (i in 1:n) {
  dens <- 1
  while (any(dens < 0.03) || any(dens > 0.97)) {
    parchoices <- seq(minval, maxval, 0.001)
    param <- sample(parchoices, 3)
    param_memory <- sample(seq(minval, maxval, 0.001), 1)  # stability [-3; 3]
    param_memory_pos <- sample(seq(0, maxval, 0.001), 1)  # pos stability [0; 3]
    param_rate <- sample(seq(5, 10, 0.001), 1)  # choose rate parameter [5; 10]
    tryCatch(
      {
        print(c(param, param_memory, param_memory_pos, param_rate))
        obj_tergm <- simul_tergm(param[1], param[2], param[3], param_memory)
        obj_tergm_pos <- simul_tergm(param[1], param[2], param[3], 
                                     param_memory_pos)
        obj_saom <- simul_saom(param[1], param[2], param[3], param_rate)
        dens <- c(sapply(obj_tergm$networks, gden)[-1], 
                  sapply(obj_saom$networks, gden)[-1], 
                  sapply(obj_tergm_pos$networks, gden)[-1])
        if (any(dens < 0.03) || any(dens > 0.97)) {
          failures <- failures + 1
          cat("Density out of bounds:", mean(dens), " (dens:", param[1], 
              ", recip:", param[2], ", transtrip:", param[3], ", memory/rate:", 
              param_memory, "/", param_memory_pos, "/", param_rate, ")\n")
        } else {
          cat("Solution:", mean(dens), " (dens:", param[1], ", recip:", 
              param[2], ", transtrip:", param[3], ", memory/rate:", 
              param_memory, "/", param_memory_pos, "/", param_rate, ")\n")
        }
      }, 
      error = function(cond) {
        failures <- failures + 1
        print(cond)
        cat("Degenerate simulations:", mean(dens), " (dens:", param[1], 
            ", recip:", param[2], ", transtrip:", param[3], ")\n")
      }
    )
  }
  sim_tergm[[i]] <- obj_tergm
  sim_tergm_pos[[i]] <- obj_tergm_pos
  sim_saom[[i]] <- obj_saom
}

save(sim_tergm, sim_saom, numnodes, timesteps, minval, maxval, n, nsim, 
    file = "sim.RData")


# ==============================================================================
# Diagnostics: print failed simulations and plot densities over time
# ==============================================================================

print(failures / (failures + n))  # percentage of failed simulations: 0.5145631

# plot density of networks over time
densities_tergm <- sapply(sim_tergm, function(x) sapply(x$networks, gden))
densities_tergm <- data.frame(Density = c(densities_tergm), 
                              Time = rep(1:(timesteps + 1), n), 
                              Iteration = sort(rep(1:n, (timesteps + 1))))
densities_tergm_pos <- sapply(sim_tergm_pos, function(x) sapply(x$networks, 
                                                                gden))
densities_tergm_pos <- data.frame(Density = c(densities_tergm_pos), 
                              Time = rep(1:(timesteps + 1), n), 
                              Iteration = sort(rep(1:n, (timesteps + 1))))
densities_saom <- sapply(sim_saom, function(x) sapply(x$networks, gden))
densities_saom <- data.frame(Density = c(densities_saom), 
                              Time = rep(1:(timesteps + 1), n), 
                              Iteration = sort(rep(1:n, (timesteps + 1))))
densities <- rbind(densities_tergm, densities_tergm_pos, densities_saom)
densities$Model <- c(rep("TERGM", n * (timesteps + 1)), 
                     rep("TERGM (positive stability)", n * (timesteps + 1)), 
                     rep("SAOM", n * (timesteps + 1)))

pdf("simulation-densities.pdf", height = 4, width = 8)
ggplot(densities, aes(y = Density, 
                      x = Time, 
                      colour = Model, 
                      iteration = Iteration)) + geom_smooth() + theme_bw()
dev.off()


# ==============================================================================
# Goodness-of-fit assessment for the TERGM process (warning: takes a while)
# ==============================================================================

tergmsim_tergm <- list()
tergmsim_saom <- list()
tergmsim_goftergm <- list()
tergmsim_gofsaom <- list()
for (j in 1:n) {
  message(paste0("Starting iteration ", j, " out of ", n, "..."))
  # estimate SAOM
  l <- lapply(sim_tergm[[j]]$networks, as.matrix)
  pred.target <- l[[length(l)]]  # last network (for out-of-sample prediction)
  training.data <- l[-length(l)]  # first five time steps (for estimation)
  last.transition <- l[(length(l) - 1):length(l)]
  siena.nets <- sienaNet(array(unlist(training.data), dim = c(numnodes, 
      numnodes, length(training.data))), allowOnly = FALSE)
  mymodel <- sienaModelCreate(cond = FALSE, useStdInits = FALSE, 
      projname = "myproject")
  mydata <- sienaDataCreate(siena.nets)
  myeff <- getEffects(mydata)
  sink("aux")  # divert screen output because includeEffects is very verbose
  myeff <- includeEffects(myeff, transTrip)
  siena.nets2 <- sienaNet(array(unlist(last.transition), dim = c(numnodes, 
      numnodes, length(last.transition))), allowOnly = FALSE)
  mydata2 <- sienaDataCreate(siena.nets2)
  myeff2 <- getEffects(mydata2)
  myeff2 <- includeEffects(myeff2, transTrip)
  sink(NULL)
  # estimate TERGM, then do GOF assessment for TERGM and SAOM
  tryCatch(
    {
      bt <- btergm(training.data ~ edges + mutual + ttriple + 
          memory("stability"), verbose = FALSE)
      tergmsim_tergm[[j]] <- bt
      mem <- training.data[[length(training.data)]]
      btgof <- gof(bt, formula = pred.target ~ edges + mutual + ttriple + 
          edgecov(mem), target = pred.target, nsim = nsim, statistics = c(esp, 
          dsp, odeg, ideg, geodesic, istar, rocpr), verbose = FALSE)
      tergmsim_goftergm[[j]] <- btgof
      sm <- siena07(mymodel, data = mydata, effects = myeff, silent = TRUE, 
          verbose = FALSE, batch = TRUE, returnDeps = TRUE)
      tergmsim_saom[[j]] <- sm
      saomgof <- gof(sm, outofsample = TRUE, sienaData = mydata2, 
          sienaEffects = myeff2, nsim = nsim, verbose = FALSE,
          statistics = c(esp, dsp, odeg, ideg, geodesic, istar, rocpr))
      tergmsim_gofsaom[[j]] <- saomgof
    }, 
    error = function(e) {
      message(paste("Iteration", j, "is dropped because the following error", 
          "occurred:", e))
    }, 
    finally = {}
  )
}

save(tergmsim_gofsaom, tergmsim_goftergm, tergmsim_tergm, tergmsim_saom, 
     file = "tergmgof.RData")

lapply(tergmsim_saom, function(x) print(x$tconv.max))


# ==============================================================================
# Goodness-of-fit assessment for the TERGM process with positive stability
# ==============================================================================

tergmsimpos_tergm <- list()
tergmsimpos_saom <- list()
tergmsimpos_goftergm <- list()
tergmsimpos_gofsaom <- list()
for (j in 1:n) {
  message(paste0("Starting iteration ", j, " out of ", n, "..."))
  # estimate SAOM
  l <- lapply(sim_tergm_pos[[j]]$networks, as.matrix)
  pred.target <- l[[length(l)]]  # last network (for out-of-sample prediction)
  training.data <- l[-length(l)]  # first five time steps (for estimation)
  last.transition <- l[(length(l) - 1):length(l)]
  siena.nets <- sienaNet(array(unlist(training.data), dim = c(numnodes, 
      numnodes, length(training.data))), allowOnly = FALSE)
  mymodel <- sienaModelCreate(cond = FALSE, useStdInits = FALSE, 
      projname = "myproject")
  mydata <- sienaDataCreate(siena.nets)
  myeff <- getEffects(mydata)
  sink("aux")  # divert screen output because includeEffects is very verbose
  myeff <- includeEffects(myeff, transTrip)
  siena.nets2 <- sienaNet(array(unlist(last.transition), dim = c(numnodes, 
      numnodes, length(last.transition))), allowOnly = FALSE)
  mydata2 <- sienaDataCreate(siena.nets2)
  myeff2 <- getEffects(mydata2)
  myeff2 <- includeEffects(myeff2, transTrip)
  sink(NULL)
  # estimate TERGM, then do GOF assessment for TERGM and SAOM
  tryCatch(
    {
      bt <- btergm(training.data ~ edges + mutual + ttriple + 
          memory("stability"), verbose = FALSE)
      tergmsimpos_tergm[[j]] <- bt
      mem <- training.data[[length(training.data)]]
      btgof <- gof(bt, formula = pred.target ~ edges + mutual + ttriple + 
          edgecov(mem), target = pred.target, nsim = nsim, statistics = c(esp, 
          dsp, odeg, ideg, geodesic, istar, rocpr), verbose = FALSE)
      tergmsimpos_goftergm[[j]] <- btgof
      sm <- siena07(mymodel, data = mydata, effects = myeff, silent = TRUE, 
          verbose = FALSE, batch = TRUE, returnDeps = TRUE)
      tergmsimpos_saom[[j]] <- sm
      saomgof <- gof(sm, outofsample = TRUE, sienaData = mydata2, 
          sienaEffects = myeff2, nsim = nsim, verbose = FALSE,
          statistics = c(esp, dsp, odeg, ideg, geodesic, istar, rocpr))
      tergmsimpos_gofsaom[[j]] <- saomgof
    }, 
    error = function(e) {
      message(paste("Iteration", j, "is dropped because the following error", 
          "occurred:", e))
    }, 
    finally = {}
  )
}

save(tergmsimpos_gofsaom, tergmsimpos_goftergm, tergmsimpos_tergm, 
     tergmsimpos_saom, file = "tergmposgof.RData")

lapply(tergmsimpos_saom, function(x) print(x$tconv.max))


# ==============================================================================
# Goodness-of-fit assessment for the SAOM process (warning: takes a while)
# ==============================================================================

saomsim_tergm <- list()
saomsim_saom <- list()
saomsim_goftergm <- list()
saomsim_gofsaom <- list()
for (j in 1:n) {
  message(paste0("Starting iteration ", j, " out of ", n, "..."))
  # estimate SAOM
  l <- lapply(sim_saom[[j]]$networks, as.matrix)
  pred.target <- l[[length(l)]]  # last network (for out-of-sample prediction)
  training.data <- l[-length(l)]  # first five time steps (for estimation)
  last.transition <- l[(length(l) - 1):length(l)]
  siena.nets <- sienaNet(array(unlist(training.data), dim = c(numnodes, 
      numnodes, length(training.data))), allowOnly = FALSE)
  mymodel <- sienaModelCreate(cond = FALSE, useStdInits = FALSE, 
      projname = "myproject")
  mydata <- sienaDataCreate(siena.nets)
  myeff <- getEffects(mydata)
  sink("aux")  # divert screen output because includeEffects is very verbose
  myeff <- includeEffects(myeff, transTrip)
  siena.nets2 <- sienaNet(array(unlist(last.transition), dim = c(numnodes, 
      numnodes, length(last.transition))), allowOnly = FALSE)
  mydata2 <- sienaDataCreate(siena.nets2)
  myeff2 <- getEffects(mydata2)
  myeff2 <- includeEffects(myeff2, transTrip)
  sink(NULL)
  # estimate TERGM, then do GOF assessment for TERGM and SAOM
  tryCatch(
    {
      bt <- btergm(training.data ~ edges + mutual + ttriple + 
          memory("stability"), verbose = FALSE)
      saomsim_tergm[[j]] <- bt
      mem <- training.data[[length(training.data)]]
      btgof <- gof(bt, formula = pred.target ~ edges + mutual + ttriple + 
          edgecov(mem), target = pred.target, nsim = nsim, statistics = c(esp, 
          dsp, odeg, ideg, geodesic, istar, rocpr), verbose = FALSE)
      saomsim_goftergm[[j]] <- btgof
      sm <- siena07(mymodel, data = mydata, effects = myeff, silent = TRUE, 
          verbose = FALSE, batch = TRUE, returnDeps = TRUE)
      saomsim_saom[[j]] <- sm
      saomgof <- gof(sm, outofsample = TRUE, sienaData = mydata2, 
          sienaEffects = myeff2, nsim = nsim, verbose = FALSE,
          statistics = c(esp, dsp, odeg, ideg, geodesic, istar, rocpr))
      saomsim_gofsaom[[j]] <- saomgof
    }, 
    error = function(e) {
      message(paste("Iteration", j, "is dropped because the following error", 
          "occurred:", e))
    }, 
    finally = {}
  )
}

save(saomsim_gofsaom, saomsim_goftergm, saomsim_tergm, saomsim_saom, 
     file = "saomgof.RData")

lapply(saomsim_saom, function(x) print(x$tconv.max))


# ==============================================================================
# Visualization of GOF results
# ==============================================================================

# extract ROC and PR values
extractrocpr <- function(goflist, name) {
  temp <- sapply(goflist, function(x) x[[7]][name])
  temp[sapply(temp, is.null)] <- NA
  temp <- unlist(temp)
  names(temp) <- NULL
  return(temp)
}

tergmsim_tergm.roc <- extractrocpr(tergmsim_goftergm, "auc.roc")
tergmsim_tergm.pr <- extractrocpr(tergmsim_goftergm, "auc.pr")
tergmsim_saom.roc <- extractrocpr(tergmsim_gofsaom, "auc.roc")
tergmsim_saom.pr <- extractrocpr(tergmsim_gofsaom, "auc.pr")

tergmsim_saom.roc.rgraph <- extractrocpr(tergmsim_gofsaom, "auc.roc.rgraph")
tergmsim_saom.pr.rgraph <- extractrocpr(tergmsim_gofsaom, "auc.pr.rgraph")
tergmsim_diff.roc <- tergmsim_tergm.roc - tergmsim_saom.roc
tergmsim_diff.pr <- tergmsim_tergm.pr - tergmsim_saom.pr

tergmsimpos_tergm.roc <- extractrocpr(tergmsimpos_goftergm, "auc.roc")
tergmsimpos_tergm.pr <- extractrocpr(tergmsimpos_goftergm, "auc.pr")
tergmsimpos_saom.roc <- extractrocpr(tergmsimpos_gofsaom, "auc.roc")
tergmsimpos_saom.pr <- extractrocpr(tergmsimpos_gofsaom, "auc.pr")

tergmsimpos_saom.roc.rgraph <- extractrocpr(tergmsimpos_gofsaom, 
                                            "auc.roc.rgraph")
tergmsimpos_saom.pr.rgraph <- extractrocpr(tergmsimpos_gofsaom, "auc.pr.rgraph")
tergmsimpos_diff.roc <- tergmsim_tergm.roc - tergmsimpos_saom.roc
tergmsimpos_diff.pr <- tergmsim_tergm.pr - tergmsimpos_saom.pr

saomsim_tergm.roc <- extractrocpr(saomsim_goftergm, "auc.roc")
saomsim_tergm.pr <- extractrocpr(saomsim_goftergm, "auc.pr")
saomsim_saom.roc <- extractrocpr(saomsim_gofsaom, "auc.roc")
saomsim_saom.pr <- extractrocpr(saomsim_gofsaom, "auc.pr")

saomsim_saom.roc.rgraph <- extractrocpr(saomsim_gofsaom, "auc.roc.rgraph")
saomsim_saom.pr.rgraph <- extractrocpr(saomsim_gofsaom, "auc.pr.rgraph")
saomsim_diff.roc <- saomsim_tergm.roc - saomsim_saom.roc
saomsim_diff.pr <- saomsim_tergm.pr - saomsim_saom.pr

t.test(tergmsim_tergm.roc, tergmsim_saom.roc)
# t = 7.2089, df = 179.65, p-value = 1.51e-11
t.test(tergmsim_tergm.pr, tergmsim_saom.pr)
# t = 4.845, df = 195.61, p-value = 2.57e-06
t.test(tergmsimpos_tergm.roc, tergmsimpos_saom.roc)
# t = 1.2959, df = 195.62, p-value = 0.1965
t.test(tergmsimpos_tergm.pr, tergmsimpos_saom.pr)
# t = 0.45891, df = 195.44, p-value = 0.6468
t.test(saomsim_tergm.roc, saomsim_saom.roc)
# t = -1.9794, df = 194.84, p-value = 0.04918
t.test(saomsim_tergm.pr, saomsim_saom.pr)
# t = -1.4037, df = 197.09, p-value = p-value = 0.162

dat <- data.frame(Value = c(tergmsim_saom.roc, 
                            tergmsim_tergm.roc, 
                            tergmsim_diff.roc, 
                            tergmsim_saom.pr, 
                            tergmsim_tergm.pr, 
                            tergmsim_diff.pr, 
                            tergmsimpos_saom.roc, 
                            tergmsimpos_tergm.roc, 
                            tergmsimpos_diff.roc, 
                            tergmsimpos_saom.pr, 
                            tergmsimpos_tergm.pr, 
                            tergmsimpos_diff.pr, 
                            saomsim_saom.roc, 
                            saomsim_tergm.roc, 
                            saomsim_diff.roc, 
                            saomsim_saom.pr, 
                            saomsim_tergm.pr, 
                            saomsim_diff.pr), 
                  Measure = c(rep("ROC", n * 3), 
                              rep("PR", n * 3), 
                              rep("ROC", n * 3), 
                              rep("PR", n * 3), 
                              rep("ROC", n * 3), 
                              rep("PR", n * 3)), 
                  Model = rep(c(rep("SAOM", n), 
                            rep("TERGM", n), 
                            rep("diff", n)), 6), 
                  Process = c(rep("TERGM DGP", n * 6), 
                              rep("TERGM (positive stability) DGP", n * 6), 
                              rep("SAOM DGP", n * 6)))

col_tergm <- "#00BA38"   # green
col_saom <- "#F8766D"    # red
col_random <- "#619CFF"  # blue

pdf("simulation-auc.pdf", width = 9, height = 4)
ggplot(dat, aes(y = Value, fill = Model, x = Measure, process = Process)) + 
    geom_hline(yintercept = 0, colour = "gray") + 
    geom_boxplot(position = position_dodge(0.85)) + 
    facet_wrap(~ Process) + 
    theme_bw() + 
    theme(legend.position = c(0.5, 0.1), 
          legend.direction = "horizontal", 
          strip.background = element_blank()) + 
    ylab("Area Under the Curve") + 
    xlab("Goodness-of-Fit Measure") + 
    scale_fill_manual(values = c(col_random, col_saom, col_tergm))
dev.off()


# ==============================================================================
# Visualization of boxplot GOF comparison results
# ==============================================================================

comparison <- list()
for (i in 1:n) {
  for (j in 1:6) {
    if (!is.null(tergmsim_goftergm[[i]][[j]])) {
      t_stats <- tergmsim_goftergm[[i]][[j]]$stats
      t_stats$names <- rownames(t_stats)
      s_stats <- tergmsim_gofsaom[[i]][[j]]$stats
      s_stats$names <- rownames(s_stats)
      m <- merge(t_stats, s_stats, by = "names", all = TRUE)
      m[is.na(m)] <- 0
      d <- abs(m$obs.x - m$median.y) - abs(m$obs.x - m$median.x)
      rn <- m$names
      sn <- rep(names(tergmsim_goftergm[[i]])[j], length(d))
      dat <- data.frame(Difference = d, 
                        Parameter = rn, 
                        Statistic = sn, 
                        Simulation = rep(i, length(d)), 
                        DGP = rep("TERGM DGP", length(d)))
      comparison[[length(comparison) + 1]] <- dat
    }
    if (!is.null(tergmsimpos_goftergm[[i]][[j]])) {
      t_stats <- tergmsimpos_goftergm[[i]][[j]]$stats
      t_stats$names <- rownames(t_stats)
      s_stats <- tergmsimpos_gofsaom[[i]][[j]]$stats
      s_stats$names <- rownames(s_stats)
      m <- merge(t_stats, s_stats, by = "names", all = TRUE)
      m[is.na(m)] <- 0
      d <- abs(m$obs.x - m$median.y) - abs(m$obs.x - m$median.x)
      rn <- m$names
      sn <- rep(names(tergmsimpos_goftergm[[i]])[j], length(d))
      dat <- data.frame(Difference = d, 
                        Parameter = rn, 
                        Statistic = sn, 
                        Simulation = rep(i, length(d)), 
                        DGP = rep("TERGM (positive stability) DGP", length(d)))
      comparison[[length(comparison) + 1]] <- dat
    }
    if (!is.null(saomsim_goftergm[[i]][[j]])) {
      t_stats <- saomsim_goftergm[[i]][[j]]$stats
      t_stats$names <- rownames(t_stats)
      s_stats <- saomsim_gofsaom[[i]][[j]]$stats
      s_stats$names <- rownames(s_stats)
      m <- merge(t_stats, s_stats, by = "names", all = TRUE)
      m[is.na(m)] <- 0
      d <- abs(m$obs.x - m$median.y) - abs(m$obs.x - m$median.x)
      rn <- m$names
      sn <- rep(names(saomsim_goftergm[[i]])[j], length(d))
      dat <- data.frame(Difference = d, 
                        Parameter = rn, 
                        Statistic = sn, 
                        Simulation = rep(i, length(d)), 
                        DGP = rep("SAOM DGP", length(d)))
      comparison[[length(comparison) + 1]] <- dat
    }
  }
}

comparison <- do.call(rbind, comparison)
comparison$Parameter <- factor(comparison$Parameter, levels = as.character(sort(
     as.numeric(as.character(unique(comparison$Parameter))))))
comparison$DGP <- factor(as.character(comparison$DGP), 
     levels = c("SAOM DGP", "TERGM DGP", "TERGM (positive stability) DGP"), 
     labels = c("SAOM DGP", "TERGM DGP", "TERGM (positive stability) DGP"))
dat_esp <- comparison[comparison$Statistic == "Edge-wise shared partners" & 
                      !comparison$Parameter %in% c("16", "17", "18"), ]
dat_dsp <- comparison[comparison$Statistic == "Dyad-wise shared partners" & 
                      !comparison$Parameter %in% c("17", "18"), ]
dat_ind <- comparison[comparison$Statistic == "Indegree" & 
                      !comparison$Parameter %in% c("19"), ]
dat_out <- comparison[comparison$Statistic == "Outdegree" & 
                      !comparison$Parameter %in% c("19"), ]
dat_geo <- comparison[comparison$Statistic == "Geodesic distances" & 
                      !comparison$Parameter %in% c("14", "15"), ]
dat_kst <- comparison[comparison$Statistic == "Incoming k-star" & 
                      !comparison$Parameter %in% c("19"), ]

comp_esp <- ggplot(dat_esp, aes(y = Difference, x = Parameter)) + 
    theme_bw() + 
    geom_hline(yintercept = 0, colour = "gray") + 
    geom_boxplot(outlier.colour = NA) + 
    coord_cartesian(ylim = c(-47, 29)) + 
    facet_wrap(~ DGP) + 
    theme(legend.position = "none", 
          axis.title.x = element_blank(), 
          axis.text.y = element_blank(), 
          strip.background = element_blank()) + 
    ylab("Edge-wise shared partners")

comp_dsp <- ggplot(dat_dsp, aes(y = Difference, x = Parameter)) + 
    theme_bw() + 
    geom_hline(yintercept = 0, colour = "gray") + 
    geom_boxplot(outlier.colour = NA) + 
    coord_cartesian(ylim = c(-103, 84)) + 
    facet_wrap(~ DGP) + 
    theme(legend.position = "none", 
          axis.title.x = element_blank(), 
          axis.text.y = element_blank(), 
          strip.background = element_blank(), 
          strip.text.x = element_blank()) + 
    ylab("Dyad-wise shared partners")

comp_ind <- ggplot(dat_ind, aes(y = Difference, x = Parameter)) + 
    theme_bw() + 
    geom_hline(yintercept = 0, colour = "gray") + 
    geom_boxplot(outlier.colour = NA) + 
    coord_cartesian(ylim = c(-3.7, 3.7)) + 
    facet_wrap(~ DGP) + 
    theme(legend.position = "none", 
          axis.title.x = element_blank(), 
          axis.text.y = element_blank(), 
          strip.background = element_blank(), 
          strip.text.x = element_blank()) + 
    ylab("Indegree")

comp_out <- ggplot(dat_out, aes(y = Difference, x = Parameter)) + 
    theme_bw() + 
    geom_hline(yintercept = 0, colour = "gray") + 
    geom_boxplot(outlier.colour = NA) + 
    coord_cartesian(ylim = c(-3.7, 3.7)) + 
    facet_wrap(~ DGP) + 
    theme(legend.position = "none", 
          axis.title.x = element_blank(), 
          axis.text.y = element_blank(), 
          strip.background = element_blank(), 
          strip.text.x = element_blank()) + 
    ylab("Outdegree")

comp_geo <- ggplot(dat_geo, aes(y = Difference, x = Parameter)) + 
    theme_bw() + 
    geom_hline(yintercept = 0, colour = "gray") + 
    geom_boxplot(outlier.colour = NA) + 
    coord_cartesian(ylim = c(-59, 35)) + 
    facet_wrap(~ DGP) + 
    theme(legend.position = "none", 
          axis.title.x = element_blank(), 
          axis.text.y = element_blank(), 
          strip.background = element_blank(), 
          strip.text.x = element_blank()) + 
    ylab("Geodesic distances")

comp_kst <- ggplot(dat_kst, aes(y = Difference, x = Parameter)) + 
    theme_bw() + 
    geom_hline(yintercept = 0, colour = "gray") + 
    geom_boxplot(outlier.colour = NA) + 
    coord_cartesian(ylim = c(-550, 18500)) + 
    facet_wrap(~ DGP) + 
    theme(legend.position = "none", 
          axis.title.x = element_blank(), 
          axis.text.y = element_blank(), 
          strip.background = element_blank(), 
          strip.text.x = element_blank()) + 
    ylab("Incoming k-stars")

pdf("simulation-endogenous.pdf", width = 10, height = 13)
grid.arrange(comp_esp, 
             comp_dsp, 
             comp_ind, 
             comp_out, 
             comp_geo, 
             comp_kst, 
             ncol = 1)
dev.off()